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Abstract 

Molecular dynamics (MD) simulations involving reactive potentials can be 
used to model material failure. The empirical potentials which are used in 
such simulations are able to adapt to the atomic environment, at the ex¬ 
pense of a signihcantly higher computational cost than non-reactive poten¬ 
tials. However, during a simulation of failure, the reactive ability is needed 
only in some limited parts of the system, where bonds break or form and 
the atomic environment changes. Therefore, simpler non-reactive potentials 
can be used in the remainder of the system, provided that such potentials 
reproduce correctly the behavior of the reactive potentials in this region, and 
that seamless coupling is ensured at the interface between the reactive and 
non-reactive regions. In this article, we propose a methodology to combine a 
reactive potential with a non-reactive approximation thereof, made of a set 
of harmonic pair and angle interactions and whose parameters are adjusted 
to predict the same energy, geometry and Hessian in the ground state of the 
potential. We present a methodology to construct the non-reactive approx¬ 
imation of the reactive potential, and a way to couple these two potentials. 
We also propose a criterion for on-the-fly substitution of the reactive poten¬ 
tial by its non-reactive approximation during a simulation. We illustrate the 
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correctness of this hybrid technique for the case of MD simulation of failure 
in two-dimensional graphene originally modeled with REBO potential. 

Keywords: molecular dynamics, mechanical failure, reactive potentials, 
coupling methodology, REBO, graphene 


1. Introduction 

The atomic interactions and structure of a material are essential in deter¬ 
mining its behavior and failure properties. Macroscopic failure mechanisms 
such as fracture or plastic deformation originate at the atomic scale, where 
events such as bond breaking and formation are the cause of cracks, vacancies 
and dislocation motions. 

Atomistic approaches, which aim at relating the macroscopic description 
of a material to the underlying microscopic motion of atoms, are appropriate 
to study the physics of mechanical failure. However, atomistic modeling of 
failure may be computationally expensive because systems of interest usually 
contain a very large number of particles (typically ~ 100, 000 atoms) and 
because events at the atomic scale are particularly complex when approaching 
failure, and thus require elaborate numerical modeling that is expensive to 
use. 

Among existing atomistic methods, quantum mechanics modeling pro¬ 
vides the most fundamental approach but requires a high computational 
cost that limits its use to small systems and small time scales (typically less 
than 1000 atoms for a maximum of a few picoseconds). In contrast, classi¬ 
cal molecular dynamics (MD) based on empirical potentials is an alternative 
technique that allows to reach system sizes and time scales relevant for the 
study of failure [1, 2, 3]. It relies on the assumption that atoms move ac¬ 
cording to the laws of classical mechanics, where the force held is derived 
from an empirical potential representing the inter-atomic interactions. Ac¬ 
cordingly, classical MD is not a ‘hrst principle’ approach, and the results of 
the simulation depend on the choice of the energy potential. 

Therefore, the key ingredient in MD is the choice of the energy potential 
to model the mutual interactions of the atoms. This choice fully determines 
the physical relevance of the simulation. However it also has a signihcant 
impact on the computational cost of the simulation. Energy potentials are 
functions of the atom (or group of atoms) positions. Classical potentials are 
valid in the vicinity of a given atomic conhguration, e.g., for solids in their 
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elastic regime close to the reference configuration. In contrast, some poten¬ 
tials, called ‘reactive’, have the ability to adapt the force held acting on an 
atom to its local environment. This extends their validity to a much wider 
range of conhgurations, e.g., mechanical failure of solids. Reactive poten¬ 
tials are usually calibrated on experiments and/or accurate quantum calcu¬ 
lations [4, 5, 6, 7, 8]. Most reactive potentials use a bond order formalism 
in which the energy potential is written as a sum of adaptive pairwise inter¬ 
actions that continuously adjust to the local environment. Since mechanical 
failure involves complex atomic rearrangements, such reactive potentials are 
often necessary in classical MD simulations of failure. In return of their 
adaptive capabilities, reactive potentials have a higher computational cost 
compared with simpler non-reactive potentials. For example, the computa¬ 
tional costs of REBO [9], Tersoff [10], AIREBO [7] and reaxFF [8] potentials, 
as implemented in LAMMPS [11], are about 18, 25, 200 and 680 times that 
of a non reactive potential made of harmonic bonds, respectively, while the 
cost of a harmonic angle is about 3.5 times that of a harmonic bond.^ 

Atomic rearrangements during failure are often very localized and the 
adaptive capability of the potential is unnecessary elsewhere. Thus, sub¬ 
stituting the reactive potential with a non-reactive analogue in all regions 
where reactivity is unneeded could signihcantly decrease the computational 
cost without affecting much the simulation results. 

Following this idea, in this article, we propose a methodology to construct 
a hybrid MD simulation in which a reactive potential is coupled to a non¬ 
reactive analogue. The reactive potential is used where the reactive ability 
is needed and the non-reactive analogue is used elsewhere, where neither the 
stress nor the temperature are high, and the material behavior is elastic. 


'^The origin of these ratios is twofold. The ratios between the performance of REBO, 
Tersoff, AIREBO and reaxFF are the results of the benchmark test provided with 
LAMMPS [11] (http://lammps. sandia.gov) for systems of 32, 000 atoms. The comparison 
of these performances with that of a harmonic bond potential comes from simulations of a 
face-centered cubic crystal with first neighbors bond interactions (12 neighbors per atom). 
No effect of system size was observed (systems of sizes between 500 and 250, 000 atoms 
were considered). Note that the performance of a harmonic bond potential directly de¬ 
pends on the number of bonds in the system: a diamond crystal with 4 first neighbors per 
atom would have a computational cost three times lower than the case considered in this 
example. Moreover, bonded interactions often include angular terms, the computational 
cost of which is about 3.5 times that of bonds because of the angle computation, which is 
more expensive than the distance computation. 
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In what follows, we will refer to the non-reactive analogue as the reduced 
potential, since it reproduces the reactive potential in a small region of the 
phase space close to a particular atomic topology and has no reactive ability. 

The success of the hybrid approach depends on: (1) the ability of the 
reduced potential to correctly reproduce the reactive one when used; (2) the 
ability to avoid or circumvent nonphysical transition effects at the interface 
between the two potentials (seamless coupling). In this work, we develop an 
approach in which the reduced potential has a harmonic formulation that 
matches exactly the second order approximation of the reactive potential 
with respect to the positions of atoms in the ground state of the material 
(i.e., when the solid is unstressed and at 0 K). Doing so, we obtain a reduced 
potential that produces a dynamics identical to that of the reactive potential 
in the vicinity of the ground state, and our implementation of a seamless 
coupling is greatly facilitated for any combination of potentials in a hybrid 
system. 

In the literature, similar hybrid techniques have been already proposed. 
For example, Buehler et ah [12] have coupled the reaxFF potential with 
the Tersoff potential to study dynamic crack propagation in a silicon single 
crystal. In contrast, we come up with a hybrid simulation technique that 
signihcantly differs from this approach (for instance we do not need a transi¬ 
tion layer since the continuity of the energy and its derivatives is ensured by 
construction). We pay special attention to the validity of the coupling from a 
methodological point of view. In particular, we investigate the compatibility 
conditions for the two potentials and the appropriate way to make the tran¬ 
sition from one potential to another. In a broader context, hybrid schemes in 
which quantum mechanics and molecular dynamics are combined have also 
been proposed [13], and they have some issues in common to MD/MD cou¬ 
pling, namely: whether or not the conservation of the energy is necessary, 
the conditions for a seamless coupling, ... 

We consider pristine two-dimensional graphene as a case study to illus¬ 
trate our methodology. Graphene is a material made of carbon atoms ar¬ 
ranged on a two-dimensional honeycomb lattice. For simplicity, without loss 
of generality, we work in two spatial dimensions, i.e., out-of-plane bending 
and vibration modes are discarded. Although this affects the physical prop¬ 
erties of graphene (in particular, phonon spectra and thermal expansion [14]), 
this is not a shortcoming for the methodology itself. 

The reactive potential we use to model pristine graphene is the 2nd- 
generation reactive empirical bond-order potential [9], REBO. This potential 
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accounts for the angular forces and gives the correct cohesive energy and 
equilibrium lattice constants of graphene, but is known to only imperfectly 
reproduce its elastic constants and phonon spectra [15]. This is a minor issue 
since we focus here on the coupling methodology. 

In what follows, we describe the molecular dynamics model of pristine 
graphene in which either the reactive potential or its reduced version is used 
(Section 2). The methodology for coupling both potentials in hybrid models 
is explained in Section 3. We next propose a criterion to switch from one po¬ 
tential to the other along the simulation, in an on-the-fly manner (Section 4). 
Results obtained by using hybrid simulations of failure are presented and dis¬ 
cussed in Section 5, followed by conclusions. 

2. Molecular model of pristine graphene 
2.1. Molecular dynamics 

We hrst describe the molecular dynamics of pristine graphene. The molec¬ 
ular system consists of a set of N carbon atoms of mass m, whose spatial 
positions {r;,r 2 , ••• ,rAr} evolve in time according to the laws of classical 
mechanics: 

yi<I<N mri = -VlE^ri,r 2 ,-,rN) (1) 

where E(^rur 2 ,-,rN) inter-atomic energy potential, a function that de¬ 

pends on the positions of the atoms. The Langevin dynamics extends these 
equations to include friction and the stochastic effect of thermal agitation: 

< I < N mri = -ViE(ri,r 2 ,-,rr,) “ imvi + ^ 27 kBTmi?(i) (2) 

where 7 is the damping constant, T is the temperature, ke is the Boltzmann’s 
constant and R(t) is a delta-correlated stationary Gaussian process with zero- 
mean, satisfying {Ra{t}) = 0 and {Ra{t)Ra{t')) = (here 5 is the Dirac’s 

delta and a represents the spatial coordinate). 

The ground state, G. S. , is the conhguration in which the total interaction 
energy reaches the absolute minimum. In the case of graphene it corresponds 
to a honeycomb lattice of parameter where d is the distance between 
the closest carbons (the value of d depends on the potential used; here, it is 
of the order of 1.40 A). 

Various energy potentials can be used to model graphene. Here, we con¬ 
sider the REBO potential [9] which is a reactive potential developed to model 
any organic matter, including graphene. We describe here the formulation 
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of Stuart et al. [7]. This potential is a reformulation of Brenner-Tersoff po¬ 
tentials [5, 6], that are in turn based on Abell’s formalism [4] in which the 
binding energy is expressed as pairwise nearest-neighbor interactions that de¬ 
pend on the local atomic environment. The total REBO interaction energy 
is given by 

N N 


pREBO _ pREBOyj 

,rjv) ~ 2-^ 2-^ ^{rj,rj,rK,■■■) 
I=l J>I 


( 3 ) 


with 


r-iREBO,/J _ 


( 4 ) 


where V-^'’ and Vl^ are, respectively, the repulsive and attractive ener¬ 
gies between atoms / and J (they only depend on the relative distance 
= \f'j — 'ri\ between atoms / and J), is the so-called bond or¬ 
der parameter, in which not only I and J are involved but also other atoms 
K. REBO can be considered as a short range potential since only a limited 
number of nearby atoms determine the interaction energy of a bond. In par¬ 
ticular, in two-dimensional graphene (in which dihedral angles do not play 
any role), only the hrst and second nearest neighbors of I and J determine 
the value of . 


2.2. Reduced potential 

Our aim is to construct a reduced potential that mimics the REBO poten¬ 
tial in the ground state of graphene (i.e., the ground states of both potentials 
share the same geometry and interaction energy) and in its vicinity. Accord¬ 
ingly, the reduced potential should reproduce as much as possible the Taylor 
expansion in powers of displacements of atoms of REBO potential around the 
ground state. The hrst order terms in the Taylor expansion are the resulting 
forces on the atoms, which vanish in the ground state. The approach we 
propose consists in additionally reproducing the second order term, which is 
given by the Hessian $ of the potential energy E: 

d‘^E 

drj^drj^ 

where I, J label atoms (i, • • • ,N) and a, (3 represent the spatial coordinates 
{x,y,z). is also called the force constant matrix, since under a harmonic 
approximation the total force in the direction a that atom I experiences when 


^ — [^lajy] — 
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atom J is displaced of ujp from its position in the ground state configuration 
is given by ^lajp ■ ujp. Many macroscopic properties of a material derive 
directly from the force constant matrix; in particular, the elasticity constants 
and the phonon dispersion curves, two critical properties for mechanics and 
fracture [16, 17]. In short, the force constant matrix characterizes the energy 
of displacements at all length scales [13]. Choosing a reduced potential that 
reproduces the Hessian ensures that the material properties essential to the 
study of failure are preserved. However, properties depending on higher order 
terms arising from REBO anharmonicity are not preserved, and therefore the 
proposed methodology is not suitable for the study of all material properties. 
We discuss these aspects in Section 5. 

To summarize, we seek a reduced potential that satishes: 


^REBO 


G.S. 


= E 


Hi 


G.S. 


= E, 


G.S. 


^^REBO ^ 

duioL Q s dujoi 

^REBO I _ I 

Ig.s. Ig.s. 


= 0 

G.S. 

= ^G.S. 


( 6 ) 


where G.S. means that the conhguration is that of the ground state (the 
same for both potentials). Many potentials E^ meet these requirements. 
The most straightforward one is: 


— Eg,s. + 2 '^2 

Instead we decided to reproduce this potential in a constructive manner, 
by parametrizing a set of harmonic springs and angles connecting pairs and 
triplets of atoms, respectively. These elements (here referred to as ‘springs’ 
and ‘angles’) are widely implemented in molecular dynamics codes (in con¬ 
trast to the potential (7)). This way of setting the potential up facilitates the 
on-the-fly construction or deconstruction and hence the potential substitu¬ 
tion, whereas in a formulation like (7) the contribution of each bond remains 
too implicit. We introduced this methodology previously in [18]. The poten¬ 
tial energy is expressed as a quadratic function of the spring lengths Vg and 
triplet angles 6a'- 


( 8 ) 

sGsprings aGangles 
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where 


-B? = (r, - (9) 

Ei = ica (». - e„J‘ 

One needs to find the topology of the springs and angles network and the 
values of the parameters Kg, Ga and a in order for Conditions (6) to 
be satisfied. 

First, we impose that all the springs and angles are at their minimum 
energies in the ground state (for any spring r^q^g = Tq;™® and for any angle 
deq,a = eusuriug that the ground state of is the ground state 

of i^nEBO view of Eqs. (8) and (9), the Hessian of the harmonic model 
can be written as 


<P 


HI 


IG.S. 


s€ springs 


Kg^! + 


a€angles 


where 


drs 

drio 


( 9rs 

\drji3 


G.S. 


and 


Ga<P. 


dOa 


A 


dria 


( Bdg \ 

V drjp ) 


( 10 ) 


(pS 


G.S. 


and only depend on the atoms that are connected by the spring s or 
the angle a and on their positions in the ground state, that is on the given 
topology of springs and angles. { is the basis of the subspace of 

Hessians that can be obtained with the given topology of springs and angles. 

By calibrating the values of Kg and Ga it is possible to approximate 
^REBOj^ g _ We follow the algebraic method proposed by Mounet [19] (used 
initially to enforce the acoustic sum rules and index symmetries on the 
Hessian of carbon allotropes obtained with density-functional perturbation 
theory). This method consists in projecting onto the subspace 

defined by We consider the Frobenius scalar product between 

two Hessian matrices (p ■ T = a jp lajp- The distance associ¬ 
ated to this scalar product is d(<P, T) = Find¬ 

ing the element ^^|g.s. of the subspace which is the closest to the Hes¬ 
sian is equivalent to projecting ^^^™°|g.s. onto the subspace. 

The basis { , ^a} can be orthonormalized with the scalar product pre¬ 

viously defined to obtain a new basis {Uo} in which, for any o and p, 
d{Uo,Up) = 6° (being 5° the Kronecker symbol). Then the best approxi¬ 
mation ^^^Ig.s. is obtained as the projection of ^^|g.s. onto the new basis. 










Figure 1: Schematic representation of the topology of springs and angles that connect an 
atom I to its first (A^l) and second {N2) nearest neighbors. 


i.e., <P^\g.s. = (^^Ig.s. • Uo ) Uo- Finally ^^|g.s. can be expressed 

in the original basis. We thus obtain the values of the stiffness constants 
Kg and Ga- Furthermore, when d (i^^|g.s., ^^|g.s.) = 0, then the prescribed 
topology of springs and angles is able to exactly reproduce the original force 
constant matrix. 

This methodology can be applied to any atomistic model, but when the 
system is a crystal (like graphene), then some symmetries exist, and as a 
consequence many of the angle and spring parameters are identical by sym¬ 
metry. Producing a Hessian consistent with these material symmetries is 
precisely the aim of the method described in [19]. 

In [18] we have followed this procedure to construct a reduced version of 
the REBO potential valid for pristine two-dimensional graphene. With four 
classes of harmonic interactions, as shown in Figure 1 and made precise in 
Table 1, our approximated model reproduces exactly the geometry, potential 
energy and Hessian of the reactive model in the ground state. Therefore, 
when atoms experience small displacements from their equilibrium positions 
in the ground state, this reduced potential accurately approximates the re¬ 
active potential. 

To check the quality of the approximation near the ground state and to 
quantify the effects of anharmonicity as the system moves away from that 
state, we compare the tensile response of bulk graphene obtained with REBO 
and with the reduced potential (see Figure 2). 

The crystalline structure was stretched in zigzag and armchair directions 
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Table 1: Parameters of the springs and angles used in the approximated model (Figure 1) 


Interaction 

Atoms connected 

Stiffness constant 

Equilibrium 
length or angle 

Spring 1 

I-Nl 

Ki 

= 39.86 eV/ A^ 

Teqp = 1.398 A 

Spring 2 

I-N2 

K2 

= 1.29 eV/ 

req .2 = 2.421 A 

Angle 1 

Nl-I-Nl 

Gi 

= 1.33 eV 

9eq,l = 27r/3 

Angle 2 

N2-I-N2 

G2 

= 6.87 eV 

9eq,2 = 7r/3 



Figure 2: Comparison of the tensile response of 2D graphene in armchair and zigzag 
directions obtained from MD simulation with REBO and the reduced potential shown in 
Figure 1 and made precise in Table 1. The vertical lines represent the strains at which 
the relative error in stress reaches 5%. The response of the reduced potential is a straight 
line in the (e, a) diagram. 


while the respective perpendicular directions were left undefornied. For each 
deformed state the energy was minimized and the stress computed. We used 
LAMMPS [11] to carry out these numerical experiments. These tests per¬ 
fectly illustrate how the approximation fails as the strain increases. However 
they just reveal the error in the forces (but not in the energy) of two of the 
many cases of lattice deformation that are present in a failure test. Therefore 
an adaptative criterion to exchange the potentials must be set up taking into 
account more situations of lattice deformation (Section 4). 

Another way of assessing the importance of anharmonicity is to use lat¬ 
tice dynamics and the so-called quasi-harmonic approximation [17] to observe 
how the phonon spectrum changes with temperature. Note that there is no 
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thermal expansion with harmonic models and that at 0 K the harmonic 
and the reduced potentials produce exactly the same phonon spectra since 
they have exactly the same Hessian. However, because of the small thermal 
expansion obtained for two-dimensional graphene (about 10“® K~^ [14]), we 
checked that only at very high temperatures the lattice dynamics differs from 
one potential to the other (see our previous work [18]). For example, the cal¬ 
culated frequency at the K point differs by less than 2 THz for a temperature 
increase of 1800 K, while the wave group velocities hardly changed. 

3. Coupling methodology 

In this section we describe the methodology to couple the reactive and 
reduced potentials in a hybrid MD simulation. 

3.1. Transition between potentials 

A possible strategy to couple the reactive and reduced potentials is to 
use a transition and a buffer zone at the interface between the reactive and 
the non-reactive regions (see e.g. Buehler et al. [12]). This is not the strat¬ 
egy we will follow but we describe it because it is a common choice. In this 
approach, both potentials are computed in the transition layer and weighted 
with a smooth interpolation function to obtain the total interaction energy. 
One could also imagine interpolating the force instead of the energy. Results 
a priori depend on the size of the transition region. The buffer zone is an 
extension of the transition zone that is needed to ensure that potentials are 
properly computed in the transition zone. Such a smooth transition ensures 
the continuity and differentiability of the hybrid potential, but raises ques¬ 
tions from a methodological point of view. In particular, any shift of energy 
between the reduced and reactive potentials (because the former is only an 
approximation of the latter) leads to artihcial forces in the transition layer 
since the weighting procedure introduces a gradient of energy orthogonal to 
the transition layer. Conversely, if forces are interpolated instead of energy, 
then there is no global underlying energy. 

Instead of a transition zone, the coupling methodology we develop here 
consists in an abrupt transition for which seamless coupling is ensured by 
considering a particular set of harmonic springs and correction forces at the 
interface between the potentials. 

Before proceeding further, we want to mention that similar problems arise 
for methods coupling reference atomistic models with coarse-grained atom- 
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istic models involving less degrees of freedom. The QuasiContinuum Method 
(QCM) is one such methods [20, 21], Several variants have been proposed. 
Some of them use a transition zone in which both models are concurrently 
taken into account through a weighting function. Some other variants are 
based on an abrupt transition between both models. The coupling can be 
either performed in terms of energy (in which case a hybrid energy is dehned 
for the whole system) or in terms of forces (in which case forces on each atom 
are dehned in a hybrid way, and there is no underlying energy for the whole 
system). We refer to [22, 23] for a review on these variants. 

We now detail our methodology. Within a hybrid model, we denote as a 
reactive bond the interaction between a pair of nearest neighbors I and J as 
formulated with the bond order formalism of Eq. (4) (the bond is precisely 
labeled by atoms I and J). Note that other atoms take part in this expression 
through the coefficient. We will dehne the reactive zone TZ as the set of 
reactive bonds and the non-reactive zone H as the set of harmonic elements. 
In the bulk of TZ and Ti zones (e.g. for the atoms R and H in Figure 3), 
interactions are modeled, respectively, with Eq. (4) or with the Eqs. (9) used 
to construct the reduced potential. 

A special treatment is needed near the interface. For instance, in Figure 3, 
the constant of the spring connecting atoms I and K is not clearly dehned 
since the bond I - J is reactive while the bond J - K is not. Indeed, with 
the topology and values of springs and angles as they are set in the previous 
section, no obvious correspondence exists between a particular reactive bond 
in Eq. (4) and a particular spring / angle in Eq. (8). Each harmonic element 
does not exclusively substitute a single reactive bond. Therefore the set of 
springs and angles near the interface has to be determined separately in order 
to ensure no redundant or missing energy between the reactive and reduced 
potentials. 

To obtain the particular values of the harmonic elements near the tran¬ 
sition interface, the methodology used in the previous section to construct 
the reduced potential can be adapted. The energy of a hybrid system can be 
written as 


77.HYB _ pHYB,REBO I pHYB,H 


( 11 ) 


with representing the total energy of the hybrid system and 

and representing the contributions of the reactive and re¬ 

duced potentials, respectively. Accordingly, the Hessian is the superposition 
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Figure 3: Schematic representation of the interface between reactive and non-reactive 
zones. Solid lines represent reactive bonds and dotted lines represent springs (for clarity, 
angles are not depicted). 


of the corresponding Hessians. We can write: = ^hyb,rebo _|_ ^hyb,h_ 

In the ground state = Esg springs + EaGangies To en¬ 

sure that the hybrid system approximates the reactive potential, the springs 
and angles constants should be chosen in order to satisfy: ^^^^|g.s. = 
(?REBO|^ g _ We do so by projecting ^ _ ^HYB,REBO|^ g 

subspace of possible ^ present case, the chosen topology of 

springs and angles allows to fully capture the reactive potential (i.e. ^^yb 
is exactly equal to This approach can be applied to any transition 

interface, but must be repeated each time the interface is modified to enlarge 
or reduce the reactive zone. In practice, one could identify and characterize 
all elementary interface geometries so as to treat any interface as the su¬ 
perposition of elementary interfaces. But doing so would be cumbersome. 
Instead, we opted for a much easier way to determine the springs and angles 
values at the interfaces, as explained hereafter in Subsection 3.2. 

3.2. Bond per bond substitution 

The bond order formalism, often used by reactive potentials including 
REBO, formulates the total energy as a sum over pairs of atoms, which we 
have referred to as reactive bonds. We take advantage of this formalism to 
propose a ‘bond per bond’ substitution technique to determine the particu¬ 
lar springs and angles constants in the vicinity of the transition interface. In 
practice, we can determine the exact set of springs and angles that approxi¬ 
mates the energy of a single REBO bond between two atoms I and J within 
a full reactive environment Indeed, we know the contribution 

single bond to the total Hessian in the ground state (because of 
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Figure 4: Schematic representation of the bond per bond substitution approach. A reactive 
bond I-J (whose energy also depends on atoms ATI to K4) can be replaced by a system 
made of seven springs and 12 angles, as made precise in Table 2. 


Table 2: Parameters of the springs and angles used in a bond per bond substitution 
approach (see Figure 4) 


Interaction Atoms connected 


Stiffness constant 


Spring 1 
Springs 2+ 
Springs 2- 

Angles 1+ 

Angles 1- 
Angles 2+ 

Angles 2- 


I-J 

I-Ki, I-K4, J-Kl, J-K2 
K1-K2, K3-K4 
J-I-Kl, J-I-K2, 

I-J-K3, I-J-KA 
K1-I-K2, K3-J-KA 
K3-I-KA, K1-J-K2 
J-K1-K2, J-K2-K1, 
I-K3-KA, I-KA-K3 


Ki = 39.86 eV/ 

K+ = 1.34 eV/ 

= -1.39 eV/ P 

G+ = 1.78 eV 

GC = -2.22 eV 
G^ = 11.51 eV 

G7 = -2.32 eV 


Equilibrium 
length or angle 
^eq,l = 1.398 A 
r+ 2 = 2.4208 A 

r-;2 = 2.421 A 

= 27r/3 

Cq.i = 27r/3 
^Jq,2 = 7^/3 

Cq,2 =7^/3 


Eqs. (3) and (5)), such that Thus, we can easily 

determine the set of strings and angles that reproduces this partial Hessian 
In the case of REBO simulation of graphene, each reactive bond can 
be replaced by a set of 7 springs and 12 angles as shown in Figure 4 and made 
precise in Table 2. We followed the same methodology explained above for 
arbitrary transition interfaces to obtain the values of the various parameters. 
In what follows, we refer to the set of springs and angles approximating a 
single reactive bond I-J as ‘harmonic bond J-J’. 

To set the hybrid system, we carry out a bond per bond substitution. 
Indeed, the springs and angles stiffness constants can be split into the contri¬ 
butions of each bond, irrespective of whether an atom is in bulk graphene or 
near the interface. Each reactive bond removed on the reactive potential side 
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adds its harmonic counterpart to the network of springs and angles. In the 
bulk of the harmonic zone, these contributions add up exactly to the values 
of springs and angles constants obtained in Section 2. In the vicinity of the 
transition interface, the total spring and angle constants differ from those val¬ 
ues. Most importantly this methodology preserves the Hessian through the 
transition interface. As an illustration, the constant of the spring connecting 
atoms J and N in Figure 3, Kj_^, can be seen as the superposition of the 
constants of the corresponding springs involved in the substitution of three 
reactive bonds {J - K, K - L and K - N), so Kj^n = K 2 = K 2 + K 2 + K 2 
(where K 2 and K 2 are given in Table 2, and K 2 in Table 1). Similarly, 
we hnd that Kk-m = -|- K 2 7 ^ K 2 since the reactive bond M-N is left 

reactive. 

The main advantage of the bond per bond substitution approach is that 
it easily adapts to any interface at almost no computational cost, thus allow¬ 
ing for fast on-the-fly modihcation of the reactive zone during a simulation. 
However, there is a remaining issue to be taken care of when the bond per 
bond substitution approach is used, namely the equilibrium of forces on the 
atoms at the transition interface. We explain the issue and the way we handle 
it in Section 3.3. 

3.3. Interface forces 

The constructed reduced potential has been set up to model systems with 
the same geometry, equilibrium energy and Hessian in the ground state as 
systems described by a full reactive potential. With both potentials, the 
system in the ground state is in equilibrium, i.e. the total forces acting on 
the atoms are equal to zero. However, the partial forces exerted by each bond 
on any atom are not necessarily zero. In contrast, by construction, all the 
forces exerted by springs and angles on the corresponding atoms are zero for 
the reduced harmonic potential; this is due to the fact that the equilibrium 
distances and angles were set a priori to the lattice distances and angles in 
the ground state. 

Therefore, both models keep atoms balanced but the underlying forces are 
different. This difference has no implication when considering bulk graphene 
entirely modeled with the reactive or the non-reactive potential, but a prob¬ 
lem arises for hybrid systems: at the interface between reactive and reduced 
potentials, atoms are not properly equilibrated since the forces caused by 
the reactive bonds near the interface are not completely canceled out by the 
set of springs and angles that make up the reduced potential. Formally, the 
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Figure 5: Illutrative example of a one-dimensional system in its ground state modeled 
with: (a) a reactive potential, (b) a reduced potential and (c) a hybrid scheme. 


Taylor expansion of the hybrid energy potential in the gronnd state exhibits 
non-zero first order terms, corresponding to the derivatives of the potential 
with respect to the interface atoms positions. This issne can be illustrated 
by the schematic example depicted in Figure 5. The same issue arises in the 
QCM method mentioned in Section 3.1. 

A strategy to solve this issue consists in adding constant balancing forces 
on the interface atoms in order to leave all the atoms balanced when the 
system is in the ground state. Equivalently, these balancing forces can be 
interpreted as new terms in the potential energy of the hybrid system, pro¬ 
portional to the displacement of the interface atoms in the direction opposite 
to the force. Since the second derivatives of these terms are zero, they do 
not affect the Hessian of the hybrid system. However, they contribute to 
the hrst order terms in the Taylor expansion of the hybrid energy potential, 
ensuring zero total values of the forces. Therefore, these balancing forces 
correct the spurious interface forces while preserving the second order ap¬ 
proximation. A similar approach has been used in some QCM approaches, 
where these additional forces are referred to as ghost forces. In practice, 
the balancing forces can be included at almost no computational cost during 
bond by bond substitution. The constant values of the forces produced by 
a single reactive bond with the geometry of the ground state are computed 
just once. Then, when at a given time in the simulation an atom I is at an 
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interface (i.e. when some but not all the bonds in which it participates are 
reactive ^), the corresponding forces, properly oriented, are added to each 
of the harmonic bonds. The forces produced by the reactive bonds together 
with these artihtial forces leave the atom force-balanced. 

4. Substitution criterion 

When the conhguration of the system of atoms is near the ground state, 
the reduced potential is an excellent approximation of the reactive potential. 
Conversely, the approximation fails when the explored configurations differ 
significantly from the ground state, because of excessive deformation or tem¬ 
perature. Therefore, we need a practical criterion quantifying the quality of 
the approximation and that can be used to trigger the substitution of the 
reactive potential by the reduced potential and vice-versa. Moreover, on- 
the-fly adaptive coupling requires a simple analytic criterion with minimum 
marginal computational cost. 

We propose to construct a criterion based on the 9 distances that fully 
characterize the environment of a reactive bond, as shown in Figure 6 (this 
environment is made of 6 atoms, i.e. 9 degrees of freedom in 2D after remov¬ 
ing the rigid body motion). These distances are computed at each timestep 
to obtain the harmonic forces. The computation of the criterion has there¬ 
fore a very limited additional computational cost. The primary objective 
of the criterion is to quantify the mismatch between the reduced potential 
and the reactive potential in such a way that values of the criterion exceed¬ 
ing a threshold value characterize appropriately excessive errors that make 
the substitution of potentials unsuitable. There is however no simple rela¬ 
tionship between the geometry and the energy or the forces. Accordingly, 
a geometry-based criterion is necessarily empirical. We evaluated different 
formulations of the criterion by assessing its relation to the error on a large 
variety of deformed configurations. In the present case, we considered the 
error in the energy only, but one could also consider the error for some other 
quantities such as the forces or the phonon frequencies of the crystal (i.e. the 
eigenvalues of the Hessian). 


^Note that in the cases of reactive potentials that involve distant neighbor interactions, 
external forces must be added not only on the atoms right at the interface between reactive 
and harmonic zones but also on each atom that is affected by the substitution of a reactive 
bond. 
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Figure 6: Schematic representation of a reactive bond I-J with the atoms involved and 
the distances that fully characterize the environment. 


Along with the expression of the criterion, one has to determine a thresh¬ 
old value. An optimal choice of the threshold is essential to ensure the quality 
of the approximation at a minimum computational cost: if the threshold is 
too small, bonds turn reactive very far from failure, and the efficiency of the 
simulation will hardly be improved; if the threshold is too large, the error 
induced by the substitution of potentials causes excessive deviations from the 
reference material behavior and spurious dynamics effects (e.g., nonphysical 
heating or cooling because of energy gaps at substitution). 

As mentioned above, a reactive bond in 2D graphene modeled with REBO 
involves 6 atoms. We constructed a criterion based on 9 distances that fully 
characterize the relative positions of those 6 atoms: the 5 nearest-neighbor 
distances (do to d^) and 4 of the second-order distances (ds to dg) (see Fig¬ 
ure 6). As a starting point, we consider as a criterion the sum over the 
9 segments of the relative differences between the distance in the current 
configuration and that in the ground state. However, we observed that the 
criterion is slightly improved (i.e. it produces a better correlation between the 
error in the energy and the criterion) when only one of the nearest-neighbor 
distances, do, is considered. Therefore the criterion we use is finally written 
as: 


with 


(7 = 1 


Cij 




G.S.I 


r 


o 
G.S. 


( 12 ) 

(13) 
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Figure 7: Deformation modes of a two-dimensional rectangle. 


We tested this criterion on a set of deformed configurations to assess 
the relation between the criterion and the error in energy for a single reac¬ 
tive bond. The set of configurations was generated as follows: we impose 
some deformations and then follow a Metropolis algorithm [24] to sample 
conhgurations according to Boltzmann’s statistics at hnite temperature, i.e., 
the temperature of interest for application of the methodology. Here, we 
consider 300 K. The specihc deformations we impose correspond to the two- 
dimensional deformation modes of a rectangle: isotropic, deviatoric, distor¬ 
tion and two bending modes (see Figure 7). The amplitude of these deforma¬ 
tions was increased linearly (in positive and negative directions for modes 1 
and 2) and, for each deformed configuration, 100 configurations were gener¬ 
ated with the Metropolis algorithm. We used a uniform jumping distribution, 
"^ 1,0 with Arj^a ~ mathcalU{—6,6) and 6 = 0.04 A, 

obtaining an average acceptance ratio of around 0.5. We expect this set 
of conhgurations to be representative of the conhgurations encountered in 
actual simulations of failure, thus reducing the bias in the testing of the 
criterion. Results are shown in Figure 8. 

With this chart, one can readily determine the appropriate criterion 
threshold corresponding to a given level of error. An ideal criterion would 
characterize the error exactly with no uncertainty. Dispersion to the right 
of the mean corresponds to conhgurations with low values of the criterion 
but large errors which is a risk for the reliability of the coupling. Dispersion 
to the left of the mean corresponds to conhgurations with large values of 
the criterion but small errors which is a potential cause of inefficiency of the 
coupling. Therefore, an appropriate criterion should keep the dispersion of 
the chart as small as possible. In fact this allows to quantify the relevance 
of a given criterion and opens the way for a more systematic investigation. 
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Figure 8: A sample of deformed configurations was generated. For each deformed config- 
uration, the values of the criterion (Eq. (12)) and of the error in the energy (-i —^—-) 
were computed. A: 2D histogram. B: Box plots of the values of the error of all the config¬ 
urations whose criterion lies below a given value Cmax- The blue whiskers represent the 
first and third quartiles, the continuous line is the median and the leftmost / rightmost 
whisker represents the 5*^ / 95*^ percentile. 


5. Graphene failure test cases and discussion 

As a test case of the proposed methodology, we performed simulation of 
failure of the 2D pristine graphene. We considered two different cases of 
fracture: tensile tests on samples with a prescribed initial crack and tensile 
tests on samples with a preexisting hole. Because of stress concentrations 
in the vicinity of the initiated crack or hole, one expects the reactive zones 
to grow from there. Stress concentrations around cracks and holes are quite 
different, thus the two test cases complement each other. 

We used an in-house implementation of the coupling methodology. The 
program LAMMPS [11] was used as a benchmark to validate our imple¬ 
mentation of classical MD and of the REBO potential. The value of the 
threshold for the criterion was set by trial and error. Periodic boundary con¬ 
ditions were imposed. For simulations with crack, the crack was initiated by 
ignoring some bonds along the zigzag direction. For simulations with hole. 
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Figure 9: Average mechanical response of the precracked system computed with the hybrid 
methodology when different thresholds are imposed to the substitution criterion. The 
dotted lines, associated to the right-hand axis, represent the fraction of harmonic atoms 
at each time. 


the hole was initiated by removing all atoms lying inside a circle. All the 
atoms were initially non-reactive, except those located where the crystalline 
structure ends (at the boundaries of the hole) or around the initially broken 
bonds. 

A Langevin thermostat was used everywhere to impose a temperature 
of 300 K. Once proper thermal equilibrium is reached, a mechanical loading 
is applied while the thermostat is kept. The systems were deformed at a 
constant strain rate until failure. The deformation was imposed in the arm¬ 
chair direction, while the system was left undeformed in the perpendicular 
direction. Atoms were not remapped during the deformation (i.e. only the 
dimensions of the simulation box were changed but the positions of atoms 
are not adjusted to the new box but are left to move according to the dy¬ 
namics). A strain rate of 0.5 % / ps was used which is small enough to 
avoid rate dependent behaviors. Therefore, we follow a quasi-static evolu¬ 
tion of the system from the ground state until failure. Post failure evolution 
is disregarded here since the choice of boundary conditions and thermostat is 
unsuitable for the study of crack propagation in brittle solids. More details 
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Figure 10: Snapshots of graphene fracture in a MD simulation with the hybrid methodol¬ 
ogy we propose (with C < 0.11) and a system with a preexisting crack. Atoms in red are 
those connected with the reactive potential. 


of these simulations are given in Appendix A. 

We compare in Figure 9 the average mechanical response (average over 
ten different realizations for each threshold) obtained when the system with 
preexisting crack is simulated with the reactive potential only and when 
it is simulated with the proposed coupling methodology. The average me¬ 
chanical response compares well over the whole loading process until failure 
appears. The mechanical behaviors start to be signihcantly different only 
with the highest threshold. Figure 10 displays some snapshots of the sys¬ 
tem during one of the simulations. Likewise, Figures 11 and 12 display the 
average mechanical responses for the system with a preexisting hole and a 
few snapshots of the corresponding molecular conhgurations, respectively. 
Again, the results of the hybrid potential accurately reproduce those of the 
reactive potential over the whole loading range. In this test differences from 
one threshold to another were less signihcant. 
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Figure 11: Average mechanical response of the system with a hole computed with the 
hybrid methodology we propose when different thresholds are imposed to the substitution 
criterion. The dotted lines, associated to the right-hand axis, represent the fraction of 
harmonic atoms at each time. 


In view of these results we can say that the coupling methodology success¬ 
fully reproduces the expected behavior of these two examples of 2D graphene 
failure. As expected, the size of the reactive zone increases as the system is 
stretched and extends in the direction of stress concentration. 

The efficiency of the coupling methodology strongly depends on the frac¬ 
tion of the system that is simulated with the reactive potential. Obviously, 
if most of the bonds are reactive, no signihcant computational improvement 
can be expected from the proposed methodology. Since the coupling method¬ 
ology adds only modest computational cost (calculation and potential substi¬ 
tution are all computation free, and the computation of interface forces only 
requires to identify interface atoms), a reasonable estimation of the ratio of 
the computational time of a hybrid simulation to that of its fully reactive 
equivalent is simply given by 

J.HYB H 

^ = (-) + ^(1 - (-)) ( 14 ) 

where and are the computational costs (seconds per atom and timestep) 
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Table 3: Values of (r) (average over the whole loading curve) resulting from the hybrid 
simulations of fracture for the different threshold values of the criterion. 

Threshold Precracked system System with a hole 


d^max 

(r) 

{r) 

0.06 

0.64 

0.55 

0.07 

0.46 

0.45 

0.11 

0.21 

0.19 


of the harmonic and reactive potentials, respectively, and (r) = ^^steps ^ 

is the average ratio of the number of reactive atoms to the total number 
of atoms N. 

The computational cost of the reduced potential is low because its im¬ 
plementation relies on lists of neighbors, i.e., the neighboring environment is 
hxed. In contrast, the reactive potential recomputes the neighboring environ¬ 
ment at each timestep. We estimated the ratio c^/c^ to be approximately 5.2 
in simulations with LAMMPS [11] on systems of two-dimensional graphene 
of 10^, 10^ and 10® atoms (no system size effect was observed). This ratio is 
rather consistent with the benchmark mentioned in the introduction^. 

The values of (r) can vary widely from one system to another. Systems in 
which stresses are highly concentrated in a small region have (r) 1 and can 

strongly benefit from the coupling ~ ^), whereas systems with very 

little stress concentration have (r) 1 so that the coupling is ineffective 

1). In the case studies presented above for validation purposes, 
stress is not so concentrated because of the small size of the system, hence 
the rather large fraction of reactive atoms (see the values of (r) reported in 
Table 3). For these test cases the computation time of the simulations may 
be reduced by a factor of three with the coupling methodology. 

The coupling methodology suffers from some limitations that we discuss 
hereafter. First, there is an inherent trade-off between efficiency and accu- 


^In introduction, we reported that REBO’s computational cost is 18 times larger than 
that of a potential made of 6 harmonic bonds per atom (face-centered cubic crystal with 
12 neighbors). In addition the computational cost of a harmonic angle is about 3.5 times 
larger than that of a harmonic bond. Therefore, for the reduced potential of 2D graphene 
which includes 4.5 harmonic bonds and 6 angles per atom, one could anticipate an efficiency 
4.5/(i8.6)+V3.5/(i8.6) = ^^^s better than REBO. 
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Figure 12: Snapshots of graphene fracture in a MD simulation with the hybrid methodol¬ 
ogy we propose (with C < 0.11) and a system with a preexisting hole. Atoms in red are 
those connected with the reactive potential. 


racy. Indeed, in the case of 2D pristine graphene modeled with REBO we 
were able to propose a reduced potential that matches exactly the Hessian 
of the reactive potential. If instead, we had selected a less complex reduced 
potential, e.g., without second neighbors springs, with a lower computational 
cost, this reduced potential would have only approximated the Hessian of the 
reactive potential. In general, reproducing exactly the Hessian of a reactive 
potential may require quite complex reduced potentials, especially if the re¬ 
active potential involves distant neighbor interactions. It is then tempting 
to consider simpler reduced potential at the expense of the ability of this 
potential to accurately reproduce the behavior of the material. For purely 
static behaviors (e.g., mechanical elasticity at low temperature), reduced 
potentials limited to close neighbors interactions are in principle sufficient. 
However, for more acute features (e.g., vibrational modes), accounting for 
distant neighbors interactions may prove to be critical. Therefore, one faces 
an inherent trade-off between accuracy and efficiency when applying the pro- 
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posed methodology. 

A second limitation of the approach is the second order approximation 
used. Indeed, various material properties, including thermal expansion and 
conductivity, non-linear mechanics, etc., depend on the anharmonicity of the 
potential, i.e., on higher order terms in the Taylor expansion. Since the re¬ 
duced potential is constructed to only reproduce the second order expansion 
of the potential, those properties are not captured by the reduced poten¬ 
tial. As a consequence, the proposed methodology is unable to estimate 
such properties. In addition, this is a source of inefficiency of the proposed 
methodology, since the substitution of potentials is valid only as long as 
anharmonic terms are small compared to harmonic terms (small strain and 
temperature). Therefore, a reactive potential with significant anharmonic¬ 
ity is substitutable only at low strain and temperature. The only way to 
circumvent this issue would be to include some anharmonicity effects within 
the reduced potential (e.g. by using non harmonic springs and angles). 

A last issue is related to on-the-fly substitution. On-the-fly substitution 
means that reactive bonds are substituted (or recreated) during the simu¬ 
lation. The differences in energy and forces between the reduced potential 
and the reactive one lead to artihcial jumps in energy and forces during that 
process. Spurious behaviors may occur. For instance, if the same reactive 
bond is substituted and recreated periodically, signihcant amounts of energy 
could be artificially pumped in or out of the system. In the present case, 
the thermostat equilibrates the system as the loading proceeds and main¬ 
tains the energy. From a practical point of view, the substitution criterion 
is critical for this issue. Substitution should occur when the error is small. 
The criterion should also avoid leaving isolated harmonic bonds inside the 
reactive regions and vice-versa. 

6. Conclusion and perspectives 

In this article, we present a methodology to speed up molecular dynam¬ 
ics simulations involving reactive potentials and, as an illustrative test case, 
we apply the proposed methodology to the simulation of failure of two di¬ 
mensional graphene. This methodology consists in setting up hybrid models 
where reactive potentials are substituted with reduced non-reactive poten¬ 
tials. The reduced potential approximating a reactive potential is constructed 
from a set of harmonic bond and angle interactions, allowing one to repro¬ 
duce the behavior of the system close to the ground state as predicted by the 
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reactive potential (same geometry, same interaction energy and same Hes¬ 
sian). The values of the parameters to use in the harmonic interactions can 
be obtained by an algebraic methodology. The reduced potential is a second 
order approximation of a reactive potential which ensures a correct lattice 
dynamics around the ground state. This approximation holds as long as 
higher order terms in the Taylor expansion of the reactive potential around 
the ground state are negligible. Non negligible contributions may be due to 
high deformation or high temperature. 

The coupling of the reactive and reduced potentials in a hybrid simulation 
is performed through an abrupt transition between reactive and non-reactive 
regions. The seamless coupling is ensured by considering a particular set of 
harmonic springs and correction forces at the interface. For potentials based 
on the bond order formalism, this particular set can be easily adapted during 
a simulation, according to the specihc reactive bonds that are being replaced 
at each time. We deal with the issue of extra forces needed to equilibrate 
the forces on the atoms at the interface between reactive and non-reactive 
regions. We solve this problem by adding constant forces, which affect neither 
the energy nor the Hessian. 

On-the-fly adaptive coupling requires a criterion to automatically switch 
from harmonic to reactive interactions during a simulation. We propose a 
criterion based on the geometrical differences between the considered con¬ 
figuration and the ground state conhguration. The precise expression of the 
criterion is calibrated by assessing its correlation with errors in energy over a 
set of plausible deformed states sampled according to Boltzmann’s statistics 
at hnite temperature. Finally we determine the criterion threshold triggering 
substitution by trial and error in hybrid MD simulations. 

We illustrate the various aspects of this methodology with a hybrid sim¬ 
ulation of fracture in graphene. The reactive potential was REBO and the 
reduced version is made of a set of springs and angles connecting each atom 
to its hrst and second order neighbors. The results of this test case satis¬ 
factorily reproduce the expected fully reactive behavior and illustrate the 
interest of the approach. The coupling methodology is efficient if only a 
small fraction of the system is simulated with the reactive potential, that is, 
if stresses and reactive processes are highly concentrated. In the example of 
the reduced potential proposed for 2D graphene the coupling can decrease 
the computational cost by a factor of 3. 

Some issues remain open to further research with the goal of improv¬ 
ing the performance of simulations without compromising the validity of the 
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results. In particular, the criterion for on-the-fly substitution can be more 
deeply studied to reduce the fraction of the system that is computed with 
reactive potentials at a given simulation time. It is also interesting to analyze 
whether or not a simpler version of reduced potential could be used far from 
the reactive zone, or how, in contrast, construct a reduced version including 
some anharmonic effects. In the end, a multi-scale hybrid approach in which 
reactive potentials and different reduced versions are coupled can make an 
important difference in the performance of molecular dynamics simulation of 
failure but it requires a careful consideration of the physical implications, es¬ 
pecially of those relating to the seamless coupling. Finally, the methodology 
presented here could be adapted to long-range reactive potentials. 
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Appendix A. Details of hybrid MD simulations 

• Features of the samples of pristine graphene: 

— Number of atoms: N = 1,800 (precracked) and N = 4,062 (hole) 

— Initial configuration: Ground state G.S. (honeycomb; a = y/Sd) 

— Distance between closest atoms in the G.S.: d = 1.3977 A 

— Potential energy per bond in the G.S.: cg.s. = = —5.2049 eV 

— Partial forces exerted by a single bond on each atom in the ground 
state (according to the schematic in Figure 6 and being x and y the 
horizontal and vertical directions and {ea;,ej^} the standard basis) 
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• General description of our hybrid MD implementation: 
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— Time integration of (2): Velocity-Verlet scheme for the Hamiltonian 
equations and exact integration of the Ornstein-Uhlenbeck process [25]: 


* ^(n+l/2) ^ ^(n) ^ 


* = -VS(,q.+i)) 



where a = exp {—^5t) and is a vector whose random compo¬ 
nents that are mutually uncorrelated, Gaussian, and have zero- 
mean and unit variance. 

— Neighbor list: Cell structures and linked cells method (see [26]) 

• Settings: 

— Timestep: 5t = 5.10“'^ ps 
— Temperature of thermostat: T = 300 K 
— Langevin damping coefficient 7 = 1.0 ps“^ 

— Number of steps for thermalization: lO'^ 

— Deformation rate 0.5 % / ps 
— Box deformation performed every step 

— Computation of the criterion / substitution performed every step 
— Neighbor list updated every 100 steps 
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